function elements = calculateLimitStage(elements,model)
ne = model.countElements();
for i=1:ne
    e = model.getElement(i-1);
    face = e.getFace();
    %check if the crack was initialized
    if elements(i).isCrack==0 && elements(i).isActivated==1
        sigma = elements(i).stress;
        %perform spectral decomposition
        oxx = sigma(1);
        oyy = sigma(2);
        oxy = sigma(3);
        t1 = (oxx + oyy) / 2;
        t2 = sqrt(((oxx - oyy)/2)^2 + oxy^2);
        o1 = t1 + t2;
        if o1 >= elements(i).ft
            disp(['warning: crack initiated at element ' num2str(i)]);
            elements(i).isCrack = 1;
            %check the positive node
            if abs(oxy) < 1e-12
                if (oxx > oyy)
                    n1 = 1;
                    n2 = 0;
                else
                    n1 = 0;
                    n2 = 1;
                end
            else
                g = sqrt((oxx - o1)^2 + oxy^2);
                n1 = -oxy / g;
                n2 = (oxx - o1) / g;
            end
            N = face.getNodes();
            cen = face.computeCenter();
            cx = cen(1);
            cy = cen(2);
            d = -n1*cx - n2*cy;
            for j = 1:3
                X = N(j).getPosition();
                f(j) = n1 * X(1) + n2 * X(2) + d;
            end
            if f(1) * f(2) > 0
                PosDomain = [3];
                NegDomain = [1 2];
            else
                if f(1) * f(3) > 0
                    PosDomain = [2];
                    NegDomain = [1 3];
                else
                    PosDomain = [1];
                    NegDomain = [2 3];
                end
            end
            %compute the crack normal direction
            posNode = N(PosDomain(1)).getPosition();
            delta(1) = posNode(1) - cx;
            delta(2) = posNode(2) - cy;
            if(n1*delta(1)+n2*delta(2) < 0)
                n1 = -n1;
                n2 = -n2;
            end
            %compute crack length
            xi1 = posNode(1);
            yi1 = posNode(2);
            negNode1 = N(NegDomain(1)).getPosition();
            negNode2 = N(NegDomain(2)).getPosition();
            xi2 = negNode1(1);
            yi2 = negNode1(2);
            xi3 = negNode2(1);
            yi3 = negNode2(2);
            a1 = yi2 - yi1;
            b1 = xi1 - xi2;
            c1 = yi1 * (xi2 - xi1) - xi1 * (yi2 - yi1);
            a2 = yi3 - yi1;
            b2 = xi1 - xi3;
            c2 = yi1 * (xi3 - xi1) - xi1 * (yi3 - yi1);
            d = -n1*cx - n2*cy;
            X1 = computeIntersection(n1, n2, d, a1, b1, c1);
            X2 = computeIntersection(n1, n2, d, a2, b2, c2);
            elements(i).Intersect1 = X1;
            elements(i).Intersect2 = X2;
            lGamma = sqrt((X2(1)-X1(1))^2 + (X2(2)-X1(2))^2);
            %compute area of element
            A = abs((xi1-xi3)*(yi2-yi1)-(xi1-xi2)*(yi3-yi1))/2;
            %compute Ge
            elements(i).Ge = -lGamma/A*[n1 0;0 n2;n2 n1];
            %compute Gc
            elements(i).crack_normal = [n1;n2];
            nb1 = (-yi2+yi3)/(xi3*(-yi1+yi2)+xi2*(yi1-yi3)+xi1*(-yi2+yi3));
            nb2 = (xi2-xi3)/(xi2*yi1 - xi3*yi1 - xi1*yi2 + xi3*yi2 + xi1*yi3 - xi2*yi3);
            elements(i).Gc = [-nb1 0;0 -nb2;-nb2 -nb1];
            %check the positive side of crack
%             if (xi1-cx)*n1 + (yi1-cy)*n2 < 0
%                 elements(i).Gc = - elements(i).Gc;
%             end
            %check the trivial relation (41)(42)
%             B = zeros(3,6);
%             B = e.computeB(face,B,[xi(1) xi(2)]);
%             xi0 = [1 1.12]';
%             dhat = zeros(6,1);
%             dhat(2*PosDomain-1) = xi0(1);
%             dhat(2*PosDomain) = xi0(2);
%             etrial = norm(B*dhat + elements(i).Gc*xi0)
%             if etrial ~= 0
%                 elements(i).Gc = -elements(i).Gc;
%             end
%             norm(B*dhat + elements(i).Gc*xi0)
            %compute Kdxi
            xi = e.getGaussPoint(0);
            w = e.getGaussWeight(0);
            jac = face.jacXAt(xi);
            detJ = abs(det(jac));
            B = zeros(3,6);
            B = e.computeB(face,B,xi);
            C = elements(i).constitutiveMatrix;
            elements(i).Kdxi = w*detJ*B'*C*elements(i).Gc;
            %compute Kxid
            elements(i).Kxid = w*detJ*elements(i).Ge'*C*B;
            elements(i).lGamma = lGamma;
            %compute Kxixi
            elements(i).Kxixi = w*detJ*elements(i).Ge'*...
                C*elements(i).Gc;
%             Ctilde = elements(i).H*[n1^2 n1*n2;...
%                 n2*n1 n2^2];
            Ctilde = elements(i).H*[n1 0;...
                n2 0];
            elements(i).Kxixi = elements(i).Kxixi+...
                lGamma*Ctilde;
            %initialize renh
            elements(i).renh = [0;0];
            %print out matrices
%             elements(i).Kdxi
%             elements(i).Kxid
%             elements(i).Kxixi
%             elements(i).Gc
%             disp(['pos domain ' num2str(N(PosDomain(1)).getId()) ...
%                 ', neg domain ' num2str(N(NegDomain(1)).getId()) ...
%                 ',' num2str(N(NegDomain(2)).getId())]);
        end
    end
end
end
